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Abstract. An effective algorithm is presented for solving 
the Beltrami equation df = ndf in a planar disk. The algo- 
rithm involves no evaluation of singular integrals. The strat- 
egy, working in concentric rings, is to construct a piecewise 
linear /i-conformal mapping and then correct the image us- 
ing a known algorithm for conformal mappings. Numerical 
examples are provided and the computational complexity is 
analyzed. 

1 Introduction 

The Beltrami equation 

df{z)/dz 
df{z)/dz 

determines a unique normalized quasiconformal self mapping / of the unit 
disk © = {2:|z|<l}in the complex plane. Here /x is a given measurable 
function in D with ||/i||oo < 1, and is called the Beltrami derivative (or 
complex dilatation) of /. We say that / is /i-conformal. 

The Beltrami equation has been studied intensively in large measure due 
to its importance in the theory of deformations of Kleinian groups and their 
applications to Teichmiiller spaces [15], [20]. Other applications of the Bel- 
trami equation are mentioned in the introduction to [5j. Some more recent 
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applications, such as mapping of the cerebral cortex use the Beltrami equa- 
tion in the spirit of its original application, dating back to Gauss, of finding 
a conformal mapping from a surface in 3-space onto a planar region; this is 
done effectively in [2J without explicit use of the Beltrami derivative /i. 

With the increasing use computer studies it has become of great interest 
to solve the Beltrami equation numerically. One method for doing this is 
suggested naturally by the classical existence proof given by Mori-Boyarskii- 
Ahlfors-Bers [I5],[20]. For this method one must evaluate singular integrals 
of the form 



(defined as Cauchy Principal Values when m = 2), and calculate sums of 
Newmann series of the form 



A related approach involving the singular integrals was developed by Daripa 
and Mashat [3]. Instead of summing the Newmann series, their method in- 
volves iteration towards a solution of a related Dirichlet problem. Their work 
incudes refinement of the technique of evaluation of the singular integrals via 
FFT, which is of interest in itself. 

G. B. Williams [32j proposed an alternative method of solving the Bel- 
trami equation, based on conformal welding of circle packings. A sequence of 
circle packings is produced whose corresponding affine mapping converge to 
the /i-conformal mapping. The emphasis in ^32j is on investigating properties 
of circle packings and on proving they can be constructed to converge to a 
solution of the Beltrami equation, and little information is given towards the 
point of view of numerical analysis. 

The present approach to the numerical solution of the Beltrami equation 
is in the spirit of 02], but with the aim of applying conformal mapping 
methods other than just circle packing. We start from the observation that 
the general solution of the Beltrami equation for constant < 1, is an 

affine linear mappping of the form 



for complex constants a, b with a ^ 0. Thus general quasiconformal mappings 
are approximated locally by such real linear affine maps. The essential idea 




m = 1,2 



f{z) = a{z + fiz) + b 
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of our method is as follows. In each element of a mesh division of a 
first attempt to construct / may be carried out by using such affine maps. 
Suppose that the piecewise linear /x-conformal mapping /o which results from 
combining these local mappings is continuous and in fact one-to-one from D 
to a domain D = /(D). Let h : D ^ 3 he a. conformal mapping, the 
existence being guaranteed by the Riemann mapping theorem. Then the 
composition f = h o Jq is a. /i-conformal self mapping of D. For numerical 
work one must usually work with Beltrami derivatives which are continuous 
or at least piecewise continuous. 

The problem of calculating a /x-conformal self-mapping of the sphere S*^ 
can also be solved numerically via the Disk Algorithm, by using domain and 
image disks covering all but an insignificant part of the Riemann sphere. Con- 
sider also the problem of constructing a conformal mapping from a surface 
in Euclidean space to D. This can be accomplished by first transforming 
it nonconformally to a disk in and then applying the Disk Algorithm to 
the induced Beltrami equation. 

In Section [2] we gather the basic facts we will need about quasiconformal 
mappings. In Section [3] we define the Disk Algorithm for solving the Bel- 
trami equation, and prove that it converges to the desired function when the 
parameters are chosen appropriately. In Section H] we discuss the choice of 
conformal mapping algorithm to be plugged into the Disk Algorithm, and 
the corresponding computational cost. Numerical examples are provided in 
Section O 

2 Preliminaries 

First we discuss affine linear quasiconformal mappings. Let fi, a, b be complex 
constants, a 7^ 0, < 1, and consider the mappings 



for 2; G C Thus L[fi] is /^-conformal and linear, while yl[a, 6] is conformal 
and affine linear. Further, the images of two points, Wi = A[a, b]{zi), i = 1,2, 
determine the coefficients a and b. All yU-conformal affine linear mappings 
are of the form A[a, b] o L[fi\, and this decomposition is unique. In particular 
we have the following. 




z + 
az + b 



(1) 
(2) 
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Proposition 1 Given Zi, Z2 distinct andwi, W2 distinct, together with \jj,\ < 
1, there is a unique fi-conformal affine linear mapping B[fi; Zi,Z2; ^1,^2] 
sending zi to wi and Z2 to W2. 

We will be interested in mappings of triangles, which can be described as 
follows. 

Proposition 2 Given Zi, Z2, z^ noncollinear and Wi, W2, W3 noncollinear, 
there is a unique affine linear mapping T[zi, Z2, z^', Wi, W2, w^] which sends Zi 
to Wi (i = 1,2,3). Its Beltrami derivative is equal to 

_ iz2 - Zi){w3 - Wi) - {z^ - Zi){w2 - Wi) 
- ^)(W3 - ^^l) - - Z^){W2 - Wi) ' 

Proof. The map T = T[zi, Z2, z^, Wi, 102,103] is well-defined because each 
of the pairs {z2 — Zi, z^ — zi), {w2 — Wi, W3 — Wi) is a basis of plane over the 
field of real numbers. To verify that the Beltrami derivative of T is (|3]), first 
consider the function 

/(C) = r[o, 1, z- 0, 1, w] (C) = -^({w- z)c + {w- z)t) 

for z,w fixed. The Beltrami derivative of / is = —{w — 'z) / {w — z) . The 
conformal linear map hi = A[\/{z2 — Z\),—Z\l{z2 — z\)\ takes 21,2:2,-23 to 
0, 1, 2 where 2 = (23 — 21)/ (22 — 21). An analogous conformal linear mapping 
/12 takes w\, W2, W3 to 0, 1, w when w is chosen appropriately. Then A is the 
composition /i^^ o f o hi and its Beltrami derivative 

as obtained from the Chain Rule, is seen to be equal to ([3]). □ 

The conformal mapping ([2]) sends triangles to similar triangles. A Bel- 
trami derivative ([3]) can be regarded as a measure of how much the triangles 
with vertices 21, 22, 23 and Wi,W2, fail to be similar. This observation will 
be used in the proof of the following technical result, which in turn will be 
used in the proof of convergence of the Disk Algorithm. 
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Lemma 2.1 Fix Zq = TqC with Tq > 0, and fix /ip, |/io| < 1- Let e > 0. 
Then there exists 5 > such that the following holds. Suppose ri > tq and 
do < d'o < di < 9[ < 02, and write 

z, = roe''', 4 = ^ie*'s fc = l,2. 

Further, suppose that \wo\ = \wi\ = \w2\ = r > and \ fi2\ < 1; <ind set 

w[ = B[fii; Zo,Zi; Wo,Wi]{z[), 

w'2 = B[fi2; Zi,Z2; Wi,W2]{z'2). 

Suppose that zi, Z2, z[, z'2 are within 6 of zq, that ^1,^2 o,re within 6 0/ /io, 
and that the ratios 

02-01 0'2 - 0[ W2-WI 

01 -0o' 01 -0o' w,-wo 

are within Sofl. Let fi be the Beltrami derivative of the affine linear mapping 
T[z[, zi, Zq; w[,Wi,Wo]. Then 

l/i - jjol < e. 

Proof. The positions of the points involved are represented schematically 
in Figure [H First suppose that zq,zi,Z2, rather than lying on the circle 
\z\ = vq, are instead collinear and further are equally spaced, and likewise for 
Wq, Wi, W2] also suppose that the segment z[, z'2 is parallel to zq, zi, Z2 and that 
1^2 ~ z[\ = \z2 — Zi\. Suppose further that fii = fi2 = /^o- One may calculate 



w'2 w[ 




Figure 1: Affine mappings of triangles 
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from Proposition [2] that the Beltrami differential of T[z[, Zi, Zq; w[,Wi,Wo] 
must be equal to fiQ. Alternatively, one may simply observe that the trian- 
gles w[,wi,wo and 102,102, wi are similar because the triangles z[,zi,zq and 
^2, 2:2, 2:1, are similar and have parallel bases. Then by elementary geometry 
the triangle wi,w'i,w'2 is similar to w[,wi,wo- Since similar to 

Zi,zi,zq, it follows again that the Beltrami derivative of the affine mapping 
is /io- 

Returning now to the hypothesis that I^^a:! = by taking 5 sufficiently 
small, which implies that ri is close to tq, one can assure that the arcs of the 
circles of radii tq, ti containing the z^, will be approximated by straight 
lines, so the Beltrami derivative /i* of the affine mapping between the middle 
triangles can be made to satisfy — /io| < e/2. Up to now we have referred 
only to equally spaced points, that is, when the ratios (jl]) are all equal to 
1. In general, since /i is a continuous function of ri, 6*1, 6*2, ^i, 6*2, /ii, and 112, 
we see that by taking 5 sufficiently small we can assure that |/x — < e/2. 
This completes the proof. □ 

We recall some well known properties of quasiconformal mappings pLj, 
PU] , [2T] . If /i , /2 are quasiconformal mappings in the same planar domain 
and have the same Beltrami derivative, then /2 = hofi where /i is a conformal 
mapping from the image of /i to the image of f2- If is measurable in 
D = {z G C: \z\ < 1} and ||/i||oo < 1, then there is a unique /x-conformal 
mapping /: D — > D satisfying the normalization 

/(0) = 0, /(1) = 1. 

In this paper all mappings will be piecewise smooth. 

3 Disk Algorithm for Solving Beltrami Equa- 
tions 

In this section we will describe the main algorithm for solving the Beltrami 
equation. By assumption a formula or program is provided for calculating 
the proposed Beltrami derivative /i(z) in all of D. 
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Figure 2: Elements of the Ring Extension. 

3.1 Ring Extension procedure 

First we describe the basic step. Let Zi, Z2, ■ ■ ■ ,Z]^ be the vertices of a pos- 
itively oriented simple closed polygon, lying within another such polygon 
whose vertices are z[^z'2^. . . ,z'j^. We suppose that the interiors of the tri- 
angles with vertices Zk, -Zfc+i, z[^-^ {k = 1, . . . , N) are disjoint (indices are 
taken modulo N). Let \n{z)\ < 1 in the doubly connected domain R between 
these two polygons. Let wi,W2, ■ ■ ■ ,wn also form a simple closed polygon in 
C. 

The Ring Extension produces w'l, W2, . . . , w'^ as follows. For each k let [ik 
be the average value of over the triangle Dk] this exists because prescribed 
Beltrami differentials are by hypothesis integrable. Define 

w'k = B[^k, Zk,Zk+i; Wk,Wk+i]{z'k), k = l,...,N, (5) 

where B is as given in Proposition [H We will say that the Ring Extension 
is applicable for the given data if the interiors of the triangles Wk, Wk+i, w^+i 
are disjoint. 

3.2 Disk algorithm 

We now define the main iteration of the algorithm for solving the Beltrami 
equation in a disk. 

The unit disk © is subdivided by a mesh with vertices Zjk, where for each 
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j, the points Zji, Zj2, . . . , Zjn form a positively oriented polygon, with 

\zjk\=rj, j = 1, . . . , M, < ri < ■ ■ • < Tm = 1- 
The arguments 

9jk = axg Zjk, k = l,...,N, 

of the vertices are determined modulo 27r. In what follows we will take these 
arguments to be uniformly spaced, though this simplification is not strictly 
necessary for the basic idea of the algorithm. 

Fix j, and assume that the images Wjk have been calculated for a /i- 
conformal mapping w = f{z) from the subdisk {\z\ < rj} to {\w\ < rj} 
normalized by /(O) = 0. We wish to calculate the images wj^i^k = 

For convenience of notation, we will suppose that 9j+i^k is close (modulo 
27r) to {9j^k + ^j,fc+i)/2 (recall Figure [1]). In all of our applications of the 
method this will be taken to be exactly so, by using the definition 

If ([6]) does not hold, then, one must modify slightly the statements and proofs 
given below. The term j/2 in the exponent of has been included in order 
to relate the indices to Figure [T] and to formula . The Disk Algorithm 
requires carrying out the following for j = 1, . . . , M. 

DISK ALGORITHM: j-th iteration. 

Step 1. Assume that the Ring Extension with Beltrami derivative /i is 
applicable to the doubly connected region Rj between Zji, Zj2, . . . , zjn 
and -2^+1,2, • • • , -2j+i,Af, mapping the inner points to the vertices 

Wji, Wj2, . . . , Wjp^. Let the resulting points be called 

According to these are specifically 

U!j+l,k = B[fJ,jk; Zj^ki Zj^k+li Wjk.,Wj^k+\\{Zj+l,k) (7) 

for A; = 1, . . . , iV. 
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Step 2. The points Wj+i^i, Wj^i^2, • • • , "W^j+i.Af, joined in order, bound a 
simply connected polygonal domain -Dj+i. Let 

hj+i-.Dj+i ^ {l^l < r^+i} 

be a conformal mapping satisfying /ij+i(0) = 0, calculated by any ap- 
propriate known method. Define 

Wj+i^k = (8) 

Step 3. Apply hj^i to each Wj/^k for j' ^ j and all k. For simplicity 
of notation we will call the images Wj>^k as well. 

Some details remain to be specified to make this general scheme precise 
and to complete the Disk Algorithm. One must decide how to start the 
algorithm by finding the first polygon wu, Wu^, . . . , Win- We have done this 
using an explicit formula [29] for the conformal mapping to D from an ellipse 
with semimajor and semiminor axes a, b (a^ — 6^ = 1) and foci ±1, 

2K 

w = \/ksio.{ sin~^M;A;^) (9) 

71 

where the Jacobi modulus k is related to the complete elliptic integral K and 
the Jacobi theta functions by 

TrK(l-m.) 

q= {a + 6)"^ = e ^(™) 




(notation from [3T]). Note that the image of the circle \z\ = ri under the 
mapping L[/i] is an ellipse with axes ri(l + |/i|), ri(l — slanted in the 
directions (1/2) arg/x, (l/2)(arg/i + 7r) respectively, modulo vr. This ellipse is 
sent by the affine map A[l/{2riy/JI), 0] to an ellipse with axes a, b. Then via 
this is transformed conformally to the unit disk and can then be rescaled 
to any desired radius, completing the first step of the Disk Algorithm. 

A more crucial decision is the choice of conformal mapping procedure in 
Step 2 of the Disk Algorithm. This will be discussed in Section HI 
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3.3 Convergence of Disk Algorithm 

Once the Disk Algorithm has been executed, the collection {wjk}j,k defines 
in a natural way a piecewise linear mapping of the polygon {zMk}k to the 
polygon {wMk}k, which sends zjk to Wjk- If desired, this mapping may be 
completed to a self mapping of D by extrapolating to the remaining points 
near the circumference but for large M this portion of the domain is in- 
significant. We now justify that this piecewise linear mapping is a valid 
approximation for the /i-conformal mapping of D. 

For purposes of normalization in the the following theorem, note that 
Zj-j/2 given by ([6]) is positive real when j is even. We write [j/2] for the 
greatest integer no greater than j/2. 

Theorem 2.1 Let \fi{z)\ < k, < 1 for all z G D, where /i is piecewise contin- 
uous in D. Let Mn, Nn oo, and for each n, let = r„o < ^ni < fn2 < ■ ■ ■ < 
i^nMn = I be chosen in such a way that supj(r„j — Tnj-i) — > as n ^ oo. 
Suppose further that the Ring Extension is always applicable so that Step 1 
can be carried out for all j . Let fn be the piecewise linear self mapping of D 
produced by the Disk Algorithm, normalized so that Wj_^/2] is positive real. 
Then {fn} converges uniformly on compact subsets of 3 to a fi-conformal 
self mapping f of D. 

Proof. Let k < k' < 1. By construction, on each outward pointing trian- 
gle of the mesh (such as z[,Zi, Zq in Figured]), /„ is an affine linear quasicon- 
formal mapping whose dilatation is no greater than k. Then by Proposi- 
tion 12.11 we have < f^' on the inward pointing triangles (such as zi, z[, z'2 
in Figure [1]) if the mesh is fine enough. Therefore fn is ^'-quasiconformal 
for sufficiently large n. By compactness of bounded families of uniformly 
quasiconformal mappings [2T], there is a limit mapping / which is also k- 
quasiconformal (it is not constant because /(O) = 0, /(I) = 1). Further, /i„ 
converges uniformly to the Beltrami differential of / almost everywhere, so 
it follows that / is /x-conformal. □ 

3.4 Condition to guarantee applicability 

We must now show how to choose radii so that the Ring Algorithm will be 
applicable in Step 1. The matter is to avoid crossing of images of adjacent 
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Figure 3: Crossing triangles in the (failed) Ring Extension. The outer polygon 
(dotted line) is not simple. 



triangles as in Figure [3l In the context of Figure [H let us say that a mapped 
triangle Wq, Wi, w[ is skewed when one of the base angles is greater than tt/2, 
thus, when the perpendicular from the apex w[ to the line containing the 
base Wo, wi of the triangle does not pass through the base. It is clear that if 
none of the triangles mapped in Step 1 is skewed, then the argument of each 
Wj^i^k will lie between those of Wjk and wj+i^k, so the Ring Extension will 
be applicable. In general, given N, we will have to take M sufficiently large, 
taking sup into account as well. This is reflected in the following result. 



Theorem 2.2 Suppose that the inequality 

^ < cos — + ^ f-TTT- sin -T7 (10 

r,- N 2|/i(^)| iV ^ ^ 

is satisfied for all z in the ring Vj < \z\ < Vj^i. Then all of the N outward- 
pointing image triangles produced in the j-th application of Step 1 of the 
algorithm are non-skewed. 



Proof. In the context of Figure [H let us suppose that Zq = e~*"z, zi = 
e''°'z, and z[ = Pz, where Vj = \z\, r^+i = Pvj, and 2a = Oi—Oq. We normalize 
the figure by subtracting the midpoint (cosq;)z of the segment from zq to Zi, 
as shown in Figure HI in which [3' = (3 — cos a. Carrying out Step 1, we take 
an average value for /i within the triangle and apply L[^]. For convenience 
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Figure 4: Proof of Theorem 

we follow this by composition with the conformal mapping A[—i/ {z — flz), 0] 
(recall the formulas ([I]), ([2])). It is thus seen that the condition for the image 
not to be skewed is 

|Rei/?'c| < sin a 



where 



z + fiz 1 + d , _ , 
a = HZ/ z. 



z — jj,z 1 — d 

Since |Rezc| = |Imc| = 2|Im(i|/|l — dp < 2|(i|/|l — dp, our condition is 
satisfied when 

^ |l-dp 



sina 2\d\ 

But a = 2tt/N and \d\ = \fi\, so we are done. □ 

The condition in this theorem is appropriate for when one is not going 
to take into account precise information regarding the distribution of \fi\ 
within ©, and in such cases it is reasonable to make the ratio (3 = rj^i/vj 
independent of j. Note that (3 = 1 + {ir/N) tana where a is the base angle 
of the isosceles mesh triangles. Since vm = 1, the radii are determined as 
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fi = 0.1 


H = 0.2 


H = 0.3 


H = 0.4 


/i = 


0.5 


32 


2.561 (5) 


1.605 (9) 


1.299 (15) 


1.156 (27) 


1.078 


(52) 


64 


1.789 (8) 


1.309 (18) 


1.155 (32) 


1.083 (58) 


1.044 


(107) 


128 


1.396 (16) 


1.156 (37) 


1.079 (70) 


1.043 (127) 


1.023 


(230) 


256 


1.198 (34) 


1.078 (80) 


1.04 (154) 


1.022 (279) 


1.012 


(504) 


512 


1.099 (71) 


1.039 (175) 


1.02 (339) 


1.011 (614) 


1.006 


(1108) 



Table 1: Ratios rj+i/rj required by the estimate (fTOl) for given values of and 
||/x||oo- The corresponding values of M prescribed by (fTT|) are shown in parentheses. 



r 



f3^~K It is also reasonable to take the central "hole" in the mesh not 
much bigger than the largest (i.e., outer) triangles, which have a base of 
l-K jN . Since r\ = [3^~'^^ , this suggests defining M at least as large as 

log/3 

if it is desired to assure fllOj) and to reasonably cover D with the mesh. (See 
Table m) Of course, larger values of M may be used if a finer mesh is desired. 



4 Conformal mapping methods and compu- 
tational cost of the Disk Algorithm 

The vast literature developed during the past century on methods for map- 
ping a simply connected domain conformally to a disk (see [S] , [H] , [15] , [30] 
and the references therein), gives an ample range of options for implementing 
Step 2 of the Disk Algorithm. 

A natural choice is to apply an osculation method [TB], for two reasons: 
the convergence is guaranteed regardless of the complexity of the original 
configuration, and most common osculation methods provide an approxima- 
tion of the mapping from non-circular domain to a disk (rather than the 
inverse direction). However, it is well known that osculation methods tend 
to converge slowly. 

The fastest available conformal mapping methods, many of which cost 
0{N log A^) computations per iteration and converge linearly or in some cases 
quadratically, arrive quickly at good approximations for the boundary maps 
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z w, and would be applicable in the Disk Algorithm because the domains 
involved to be mapped are are close to circular. However, such methods 
provide directly only the boundary values of the conformal mapping, so one 
must then use integration to calculate each inner value as a boundary integral. 
Further, it is the inverse mapping w z which is needed for the Disk 
Algorithm. Thus additional programming would be required. Since N may 
be potentially large, we have ruled out Schwarz-Christoffel methods [3], [TU]. 

In the examples computed in Section [5] we have used the Zipper Algo- 
rithm [22] ■ Briefly, this algorithm approximates a conformal mapping from 
a the interior of an A^-point Jordan curve polygon to the unit disk, as the 
composition of A^ partial zipper mappings, each of which is in turn a compo- 
sition of a few Mobius transformations and square roots. To find each partial 
zipper mapping requires requires 0{N) arithmetic operations, thus giving an 
operation count of 0(A^^) for calculating the zipper mapping. The Zipper 
Algorithm is similar to the classical osculation mappings in the sense that 
it provides an explicit formula for the approximated conformal mapping in 
the interior as well as on the boundary of the domain. It is different in that 
for fixed A^, it provides a single mapping rather than an infinite convergent 
sequence of mappings. 

The calculation of the computational cost of the Disk Algorithm with 
zipper mappings is as follows. In the jth iteration. Step 1 requires 0{N) 
arithmetic operations. Step 2 requires 0{N'^) operations to find the con- 
formal mapping, and in the process produces the A^ image points of the 
boundary. Over M values of j, Steps 1 and 2 cost 0{MN'^). In Step 3, the 
mapping, which costs 0{N) operations to map a single point, must be ap- 
plied to j collections of A^ points each, making 0{jN'^) operations. Summing 
these costs over 1 < j < M gives 0{M'^N'^) as the total operation count of 
Step 3. 

The preceding analysis is for the cost applied to the complete triangular 
mesh in ©. However, there are situations when one may only be interested in 
calculating f{z) for a reduced subset of interior points; then it is not necessary 
to apply the conformal mappings to the mesh points other than these. Then 
the cost of Step 3 may be much lower. For example, with a submesh of 
M' concentric rings of A^' points (M' < M, A^' < A^), the calculation count 
is reduced to 0{MN^) + 0{M''^N'N). As an extreme case, when one only 
wants to find the boundary map it is not necessary to keep any of the 
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inner image data; in other words, Step 3 may be suppressed entirely, and the 
operation count is then 0{MN'^). Examples of proper submeshes {M',N') 
will be given in Section [51 

Corresponding operation counts may be obtained similarly if other con- 
formal mapping algorithms are used in Step 2 of the Disk Algorithm. 

5 Numerical Results 

The Disk Algorithm was programmed in Mathematica (Wolfram), on a stan- 
dard laptop computer of approximately IGH, with the conformal mapping 
step being carried out by calling Fortran routines obtained from the website 



The following examples illustrate the results of the Disk Algorithm for a 
variety of prescribed Beltrami derivatives fi. 

(1) Constant /i. We take ^{z) = 0.4. The Disk Algorithm was applied for 
a mesh of equilateral triangles with {M,N) = (64, 128). Figure [5] shows the 
image of a subdomain with (M, A^) = (32, 32). 

Note that the image triangles in Figure O are not similar, even though /j, 
is constant, since the similarity class of an image triangle depends upon both 
the value of /i and the slope of the base of the domain triangle. 

Note that constant fi doesn't imply that all the triangles are skewed the 
same way, it depends on the position of the domain triangle too. 

(2) Radial quasiconformal mapping. Let ip: [0, 1] — > [0, 1] be an increasing 
diffeomorphism of the unit interval. Then the radially symmetric function 



of m 



,t arg 2 



(12) 



has Beltrami derivative equal to 



/i(z) 



z\{y,'{z)/^{z)) + l 



(13) 



when 2; 7^ 0. As an illustration we will take 



Vir) = 1 + -sm'^7r{r - -) 
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Figure 5: Image for constant Beltrami derivative fj, = 0.4. The original triangula- 
tion of unit disk is shown at left. 




as in Figure O 

The resulting Beltrami derivative satisfies \\fi\\ = 0.33 approximately 
(check). The Disk Algorithm was applied with {M,N) = (256,128), with 
only the real values being conserved in Step 3. These real values were com- 
pared with the true values if{r) and found to have an error no greater than 
0.0025. The image region is depicted in Figure [61 

(3) Sectorial quasiconformal mapping. Let ip: [0, 2tt] [0, 2tt] be an increas- 
ing diffeomorphism. Write -^(e*^) = e^^^^\ Then the sectorially symmetric 
function 

/(^) = klV^(^) (14) 

has Beltrami derivative equal to 



when z ^ 0. As an example we will take 



|, 0<^<7r, 
1 + ^^, n<9<2n. 



The Disk Algorithm was applied with {M,N) = (256,128), conserving no 
interior values in Step 3. The final boundary values were compared with 
the true values ip{6) and found to have an error no greater than 0.055. The 
image region is depicted in Figure [71 

(4) Exterior mappings. In Daripa [5J, quasiconformal mappings from D to 
the exterior of an ellipse (the origin being sent to oo) are calculated with the 
following two sample Beltrami derivatives, 

^,{Z) = |;,|2e0-65(-^-2.0)^ 

fJ'2{z) = -1^1^ sin(5Re z). 

Daripa uses M radii equally spaced in [0, 1], in contrast to the exponential 
spacing we have been using. The exterior mapping results can be related to 
those of the Disk Algorithm by use of the formula 

^ (l + a)-(l-a)z^ 
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Figure 7: Angular function ip of example (3 (left), together with image domain 
under sectorial mapping. 

which transforms D conformally to the exterior of an elhpse with aspect ratio 
a. Composition of h following the quasiconformal self mapping of D provides 
a mapping to the exterior of the ellipse with the same Beltrami derivative. 

To match the examples in [5J, a = 0.6 is specified. (However, the inner 
ellipses in |5] appear to have aspect ratios of approximately 0.47; axes are 
not drawn.) Our results are depicted in Figure El The images we obtain with 
(M, A^) = (64, 64) and (M, A^) = (256, 256) turn out to be indistinguishable. 
These images appear a bit different from the figures shown in [5], and it is 
not clear that the discrepancies can be accounted for simply by a different 
selection of level curves. 

Computation times are reported in |5j for = 64 as approximately 8.5 
seconds of CPU on a MIPS computer described as "approximately 15 times 
slower than the CRAY-YMP at Texas A & M University." Our laptop CPU 
times, using Mathematica/Fortran, were approximately 85 seconds. 

(5) Quasiconformal deformation of Fuchsian groups. In this example ix 
is taken to be a quadratic differential for a Fuchsian group G (see [20] 
for definitions). The elliptic modular group is the collection of conformal 
transformations of the upper half plane generated by the two mappings 
T ^ T + 1, T ^ This Fuchsian group has a standard fundamen- 
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tal domain {r : Imr > 0, iRerj < 1/2, |r| > 1}. The Klein modular 
function J satisfies the invariance relations J(t) = J(r + 1) = J(1/t), from 
which it follow that 

J'(7(r))7'(r) = J'(r) 

for every element r in the elliptic modular group. The upper half plane is 
transformed into the unit disk D via the elementary conformal mapping 



Z = h{T) 



1-T 



T 






Figure 8: Self mapping of unit disk with Beltrami derivative fii (left), and image of 
grid near unit circle under exterior conformal mapping for ellipse. Above: Beltrami 
derivative /xi; below: Beltrami derivative /X2. 
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Let G denote the corresponding group of self mappings of D. Let ip{z) = Jij). 
Then for < a < 1 the function 



is a Beltrami differential for the group G; that is, n{'^{z))'^{z) j-^^z) = n{z) 
for every 7 in G. This implies that any /i-conformal self mapping /: D ^ D 
satisfies f{j{z)) = 'y^{f{z)) where 7 i— > 7^^ is an isomorphism from G to 
another Fuchsian group G^. 

We have taken a = 0.4. Figure [9] shows part of the tesselation of D by 
fundamental domains of G and their images, which are fundamental domains 
of Gf^. Some inaccuracies are clearly visible inasmuch as the image curves 
must be hyperbolic geodesies. 



Figure 9: Self mapping of unit disk with Beltrami derivative defined by the Klein 
modular function. Computations were made with {M,N) = (256,512), {M',N') = 
(128,64). The standard fundamental domain in the upper half plane has been 
mapped to the unit disk (lighter contour); the image under the quasiconformal 
mapping is superimposed. 



fj.{z) = a 
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6 Discussion and Conclusions 



It is stated in Daripa [5] that prior to that article there were no constructive 
methods pubhshed for solving the Beltrami equation numerically. We discuss 
here some aspects of [5j in relation to the Disk Algorithm. 

As we mentioned in the Introduction, the method of is based on evalu- 
ation of singular integrals. That method is presented in the context of finding 
a /i-conformal mapping to a prescribed star-shaped domain, and is in some 
ways reminiscent of the classical method of Theodorsen [16] for conformal 
mappings. Convergence proofs are not supplied. 

Daripa's main algorithm requires evaluation of the d/dz derivatives which 
appear in the singular integrals. A variant is also proposed which does not 
require these derivatives; however this is not used in the numerical exam- 
ples provided. The operation count of one iteration of Daripa's method is 
0{MN log A^). This should be multiplied by the average number of iterations 
required, which depends on how refined the mesh is and how much accuracy 
is desired. 

The examples we have given in Section [5] show that our method begins 
to have difficulties when \n\ is large as 0.5. Likewise, in the examples which 
we have taken from |5], ||/i||oo is approximately 0.5, but it should be noted 
that \fi{z)\ is in fact bounded by 0.12 for \z\ < 0.5, and by 0.05 for \z\ < 0.3. 
In fact, an important limitation stated in [5J is that the Beltrami derivative 
/i must be Holder continuous. Further, it is recommended that /i vanish at 
least as fast as at the origin. The Disk Algorithm, in contrast, is not 
subject to any continuity requirement on fi. 

We believe that the Disk Algorithm for solving the Beltrami equation is 
conceptually simpler than other methods which have been presented, and is 
easier to implement. 
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